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Biroli et al.'s extension of the standard mode-coupling theory to inhomogeneous equilibrium states 
[Phys. Rev. Lett. 97, 195701 (2006)] allowed them to identify a characteristic length scale that 
diverges upon approaching the mode-coupling transition. We present a numerical investigation of 
this length scale. To this end we derive and numerically solve equations of motion for coefficients in 
the small q expansion of the dynamic susceptibility Xq(k; t) that describes the change of the system's 
dynamics due to an external inhomogeneous potential. We study the dependence of the characteristic 
length scale on time, wave-vector, and on the distance from the mode-coupling transition. We verify 
scaling predictions of Biroli et al. In addition, we find that the numerical value of the diverging 
length scale qualitatively agrees with lengths obtained from four-point correlation functions. We 
show that the diverging length scale has very weak k dependence, which contrasts with very strong 
k dependence of the q — > limit of the susceptibility, Xq=o(k;t). Finally, we compare the diverging 
length obtained from the small q expansion to that resulting from an isotropic approximation applied 
to the equation of motion for the dynamic susceptibility Xq(k;t). 

PACS numbers: 



I. INTRODUCTION 

As a liquid is cooled, its dynamics not only gets slower 
but also becomes increasingly heterogeneous [TJ [5J |3]. 
Moreover, the characteristic size of regions with dynam- 
ics both significantly faster and significantly slower than 
the average dynamics grows upon cooling. This obser- 
vation has led to the definition of a dynamic correlation 
length that measures the size of these so-called dynamic 
heterogeneneities. The dynamic correlation length was 
defined in terms of a four-point correlation function [4 
or a corresponding four-point structure factor [5J [HI EJ [S] ■ 

While various four-point functions can readily be ob- 
tained from simulations (albeit they typically require 
more computational effort than the familiar two-point 
functions), they are difficult to access experimentally. To 
the best of our knowledge, four-point functions have been 
obtained directly only from experiments on granular sys- 
tems [5J HO]. In a remarkable development, Berthier et 
al. [IT] showed that derivatives of standard two-point 
functions with respect to thermodynamic variables like, 
e.g., density or temperature, could be related to integrals 
of three-point correlation functions. This opened a door 
to experimental investigations of the overall degree (the 
strength) of dynamic heterogeneity upon approaching the 
glass transition jUj. It should be emphasized, however, 
that Berthier et al. could obtain the dynamic correlation 
length characterizing the spatial extent of dynamic het- 
erogeneneities only by using additional assumptions that 
related the integrals of various three-point functions to 
characteristic length scales exhibited by these functions. 

While theoretical understanding of the derivatives of 
two-point functions with respect to thermodynamic vari- 
ables is limited, Biroli et al. [13 showed that the mode- 
coupling theory could be used to analyze a closely related 
quantity, a three-point susceptibility x q (k;f), which de- 



scribes the change of the intermediate scattering function 
due to an inhomogeneous external potential. The advan- 
tage of this approach is that it allows one to evaluate a 
characteristic length scale which, up to that time, had 
remained hidden within the well-known mode-coupling 
framework. This was possible due to the fact that Biroli 
et al. considered a non-uniform external perturbation 
rather than a uniform change of density or temperature. 

To analyze the three-point susceptibility Biroli et al. 
extended the standard mode-coupling theory to describe 
the time evolution of the intermediate scattering func- 
tion of a system under the influence of an inhomoge- 
neous external potential. They defined the three-point 
susceptibility x q (k; t) as a derivative of the intermediate 
scattering function with respect to a Fourier component 
of the external potential, U(q). Biroli et al. showed that 
upon approaching the mode-coupling transition both the 
q — > limit of the three-point susceptibility and a charac- 
teristic length defined through the small g-dependence of 
X q (k; t) diverge. Moreover, they derived scaling predic- 
tions for the time-dependence of the characteristic length, 
and they found that this length grows as W 2 in the early 
P regime (here a is the mode-coupling exponent describ- 
ing approach of the intermediate scattering function to 
its plateau value) and then saturates in the late fi and 
a regimes. This was contrasted with the time depen- 
dence of the q — > limit of the three-point susceptibility 
which grows as t a and t b in the early and late (3 regimes, 
respectively (here b is the so-called von Schweidler expo- 
nent of the mode-coupling theory describing departure 
of the intermediate scattering function from its plateau 
value), peaks around the a relaxation time, and then de- 
cays to zero. The strikingly different time-dependence of 
the characteristic length and the q — > limit of the three- 
point susceptibility was interpreted as an indication of 
changing fractal dimension of dynamic heterogeneities. 
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Biroli et al. derived only scaling predictions for the 
three-point susceptibility x q (k;i). They verified some 
of these predictions through a numerical analysis of a 
schematic model that completely disregards k depen- 
dence. Here we present a numerical investigation of the 
small q behavior of the three-point susceptibility and the 
associated characteristic length. We focus on the time 
and k dependence of the length, and on its dependence 
on the distance to the mode-coupling transition. 

We start with a definition of the three-point suscepti- 
bility Xq(k; t) m Sec.jn] Next, we postulate an expansion 
of the three-point susceptibility Xq(k; m powers of q 
and derive equations of motion for the first few coeffi- 
cients in this expansion (see Sec. III). We also present 



an alternative approach to a numerical evaluation of the 
characteristic length which is based on an isotropic ap- 
proximation to the equation of motion (see Sec. IV I. 
Next, in Sees. |V| and [VTT| we describe the results of the 
numerical calculations based on the small q expansion 
and the isotropic approximation, respectively. We finish 



with a summary and conclusions in Sec. VIII 



II. THREE-POINT SUSCEPTIBILITY X q (k;t) 

To obtain the equation of motion for the three-point 
susceptibility Biroli et al. |13j considered a Newtonian 
fluid subject to a periodic in space external potential, de- 
rived mode-coupling equation of motion for the interme- 
diate scattering function of this inhomogeneous system, 
and then differentiated this equation with respect to the 
external potential. Subsequently, one of us has derived 
the equation of motion for the same three-point suscep- 
tibility for a Brownian system [2]. Not surprisingly, the 
overdamped limit of the equation of motion derived by 
Biroli et al. coincides with the equation of motion de- 
rived starting directly from Brownian dynamics. In this 
work we will use the latter equation. 

For a system subject to a non-uniform external poten- 
tial the intermediate scattering function is not diagonal 
in the wave- vector, 

F(k 1 ,k 2 ;t) = -i( /0 (k 1 ;tM-k 2 )) £/ . (1) 

Here p(ki; t) is the Fourier transform of the microscopic 
density, 



-iki-rj(t) 



(2) 



with Tj(t) being the position of the jth particle at time t. 
Furthermore, in Eq. |lj p(ki) = p{k-i',t = 0) and 
denotes the equilibrium average for a system subject to 
a static non-uniform external potential U. 

The three-point susceptibility Xq(k; t) is defined 
through an expansion of the intermediate scattering func- 
tion in powers of a harmonic external potential C/ q = 



F(k, k i; t) = F(k; t)5 kM + Xq (k; t) (-f3U ) <5 k+ q, kl + .... 

(3) 

For a system with Brownian dynamics, the equation of 
motion for the three-point susceptibility x q (k; t) has the 
following form: 

<9*X q (k; t) + ^yXq( k ; *) 

+ [ dt'M hl (k;t-t')d t , X ^;t') 
Jo 

+ J dt'MX(k;t-t')d t ,F(\k + q\;t') 

= 5 q (k;t) (4) 

In Eq. p} Dq is the diffusion coefficient of an isolated par- 
ticle, S[k) denotes the static structure factor, M m (k;t) 
is the irreducible memory function of mode-coupling the- 
ory, 



M' m {k;t) = 
uDq f dki 



(5) 

[tfe(ki,k-ki)] 2 F(*i;t)F(|k-ki|;t), 



>3«k(ki,k- ki) 



(6) 



2 J (2tt) 2 ' 
and M x (q, k; t) is defined as follows, 

M*(k-t)- nD ° k ( ^ 

xxq(ki;*)f(|k - ki|;i)«k+q(ki + q,k - ki). 

In Eqs. gg tfctkx.ka) - k • (kicfa) + k 2 c(A: 2 )) with 
k = k/fe, and n is the density. The source term in Eq. 
Q, S q (k;i), is given by 

5 q (k;t) = 



(7) 



D k 2 S(q) 1 



k - (k 



-qi; 



(k+q) 



q|;*) 



^F(|k+q|;t')- 



+S(q) / dt'M m {k;t~t')- 

Finally, the initial condition for the three-point suscep- 
tibility is x q (k; t = Q) = S(k)S(q)S(\k + q|). This form 
of the initial condition is obtained by applying a con- 
volution approximation to the exact expression for the 
initial condition, which involves a three-particle correla- 
tion function. One should note that the same convolu- 
tion approximation is used in the derivation of the mode- 
coupling equation of motion (both in a uniform and a 
non-uniform equilibrium state). 

It should be emphasized at this point that solving 
Eq. Q numerically is considerably more involved than 
solving the uniform equilibrium mode-coupling equations 
[HI [T6j [17] , and, to the best of our knowledge, has never 
been attempted. The reason is that while q is a param- 
eter in Eq. Q, non-zero value of q breaks rotational 
symmetry. Thus, Xq(k; t) depends not only on k = |k| 
and q = |q|, but also on the angle between k and q. Most 
importantly, the latter angle is an independent variable, 
rather than a parameter in Eq. Q. 
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III. EXPANSION OF Xq (k;t) 



obtain the following equation, 



A. Preliminaries 

In this work we focus on the characteristic length de- 
fined through the small g-dependence of the three-point 
susceptibility Xq(k;i). Thus, to calculate this length we 
only need to obtain the small q behavior of the suscep- 
tibility. We postulate the following expansion of Xq(k; t) 
in powers of q, 



Xq(M 



X (0) (M)+X>a 



SXq(k;i) 



dq a 



q=0 



q a 9/3 



3 2 Xq(M 



q=0 



+ E 2 L dqadqp 

X ^(k;t)+k- qx {1 Hk;t) 
Wx {2 \k;t) 

3fk.q) 2 -q 2 ) X ( 1 2) (fc;i) + ... (8) 



where the second equality follows from symmetry consid- 
erations and quantities x^\ X^ 2 \ an d x[? are defined as 
follows: 



x (1) (M) 
x (2) (M) 



x<Z\k;t) 



5> 



% q (k;t) 



dq a 



1 v 9 2 x q (k;t) 



dq a dq 



q=0 



q=0 



\Y1 [ kak p - 



a/3 



5 2 x q (k;t) 



i a q/3 



(9) 
(10) 



(11) 

q=0 



We show in Sec. [vjthat the first order term, X^{k; t), 
does not lead to a diverging characteristic length scale. 
Furthermore, we show that the second order term origi- 
nating from the trace of the second derivative of Xq(k; t), 
X^'(k', t)t leads to a diverging characteristic length scale. 
Finally, it can be shown that the second order term orig- 
inating from the symmetric traceless part of the second 

(2) 

derivative of Xq(k; t) , \ti (k; i) , does not lead to a diverg- 
ing characteristic length scale thus we omit the equation 
of motion for x\i \k;t) for sake of space a clarity. 



d tX ^(k;t) + D ---x m (k;t) 

+ / ^'Af irr (fc;i-i')9 t 'X (0) (fc;*') 



+ / dt'Mg(k;t - t')d t ,F(k;t') 



nD k 2 S(0)c(k)F(k; t) 

+5(0) / dt'M ilT (k;t-t')d t >F(k;t') (12) 
Jo 



where 



M*(k;t) = nD J ^ K(k X) k - k x )] 2 

x X < >(fc i; t)F(|k-ki|;t). (13) 

Furthermore, taking q — > limit of the initial condition 
X q (k;t = 0) = S(k)S(q)S(\k + q|) we obtain the initial 
condition for x (0) (M), X (0) (M = 0) = S(0)S(k) 2 . 

The zeroth order coefficient, X (k\t), satisfies essen- 
tially the same equation of motion as three point suscep- 
tibility Xn{k',t), which is defined as the density deriva- 
tive of the intermediate scattering function [18 . This 
was expected: in the long wavelength, q — > 0, limit the 
derivative with respect to the external potential differs 
from the derivative with respect to the density by a ther- 
modynamic factor proportional to (dn/d(3(i) T . 

The equation of motion ( [12] ) can be solved using the 
static structure factor S(k) and the dynamic scattering 
function F(k;t) as input. We calculate F(k;t) using the 
mode-coupling theory; the equation of motion for F{k; t) 
reads, 



d t F(k-t) 



S{k) 



F(k;t) 



+ / dt'M 1TI (k;t-t')d t >F(k;t') =0,(14) 
Jo 



where M 1 " is the irreducible memory function given by 
Eq. {jj). 



B. Zeroth order coefficient x (k;t) 

To get the equation of motion for X^°\k; t) we need to 
take q — > limit in all terms in Eq. Q. In this way we 



C. First order coefficient \^(k;t) 

To get the equation of motion for X ^ {k\ t) we need to 
expand all the terms in Eq. Q in powers of q and then 
to collect terms linear in q. After exploiting rotational 
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symmetry we get the following equation of motion 
+ f dt'M iTI {k;t-t')dt<x {1) {k;t') 



dt'M?(k;t-t')d t <F{k;t') 



+ / dt'M*(k-t-t')d v d k F(k-t') 
Jo 

D S(0)F(k;t)k 2 [ 1 dS(k) 1 
S(k) [lS(kf)~dk k 

+nD k 2 S(0)c(k)d k F(k;t) 



+5(0) / dt'M { "{k;t-t')d t .d k F{k;t') 



and 



(15) 



5(0) 



dt'M i "(k; t — t')d t 'F(k; if), 



where 
M*(k;t) = 

nD Q 



dki 



[Vk(k 1} k-ki)] 



■ k-ki 

kk± 



(27T) 3 1 

xFdk-kiljt)^^;*) 

+«A> I ^v w (k 1 ,-k--k 1 )F(\^--k 1 \;t) X {0) (k 1 ;t) 



x ^ c(fci) + 



[k-ki] 2 rfc(fci) ^(ki,k-ki) 



(16) 



Furthermore, we obtain the following expression for 
the initial condition for x (k[t), x^{k;t — 0) = 
S(k)S(0)dS(k)/dk. 

To solve Eq. ( 15 1 we need the equation of motion for 



the partial derivative of the intermediate scattering func- 
tion with respect to the wave-vector, d k F(k;t). This 
equation can be derived from the mode-coupling equa- 
tion ((Til: 



d t d k F{k;t) + ^^d k F{k;t) 



Dpk 
S(k) 



S(k) 

2 ~W) dkSik) 



F(k;t) 



dt'M [ "(k;t-t')d t >d k F(k;t') 



Ml l {k;t-t')dt'F{k;t') 



+ / dt'M^(k; t - t')d t >F(k; if) = 



where 
M\ x = 

nD 



(17) 
(18) 



dki . ,, ., f , , [k • p] 2 dc(p) 1 

Wf [Mkl ' p)] \ c{p) + -^r^r\ 



M$\k;t) = 



(19) 



In Eqs. l[TJ[T9]l p = k - ki and p = |k - ki 



D. Second order coefficient \^\k\t) 



For symmetry reasons, there are two linearly inde- 
pendent second order coefficients, the coefficient propor- 
tional to the trace of the matrix of second derivatives of 
X q (k; t), X^(fc; t)i an d the coefficient proportional to the 
symmetric traceless part of the matrix of second deriva- 

(2) 

fives of Xq(k; t), x\i \k\ t). It can be shown that only the 
former coefficient leads to a characteristic length that di- 
verges upon approaching the mode-coupling transition. 
Therefore, since the focus of this work is the diverging 
characteristic length, we will only give the equation of 
motion for x^(k;t). By expanding equation of motion 
([I 1 in powers of q, collecting the second order terms and 
taking a trace of the corresponding tensorial equation of 
motion we can derive the following equation of motion 
for x (2) (k;t): 



d t x^(k;t)+D ^--x i2 Hk;t) 
S(k) 

+ / d^ irr (fc;i-i')<9*'X (2) (M') 



+ / dt'M*(k;t-t')dt>F(k;t') 



dt' 'M*{k;t - lf)dt>d£F(k;lf) 
t 

dt'M${k;t-t')d v d k F(k;t') 



6 J 



3k J 



3 7o 



dt' M?(k;t - t')d t 'd k F(k;t') 
dt'M${k;t-t')d t ,F{k]t') 



xF(k i; t)F(p;t), 



+S&(k;t), 



(20) 
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where 
S^(k;t) 



(21) 



theory equation of motion, Eq. (14 1: 
d t d 2 F(k;t) 



~d 2 S(q) 




nD k 2 c(k)F(k;t) 


dq 2 


q=0_ 


2 



2D AD k dS(k) 



D k 2 S(Q) 
35(fc) 2 

D k 2 S(0) 
65(fc) 2 

d 2 S(k) 
dk 2 

1 
3 

1 5(0) 
6 k Jo 

rt 



d k F{k;t) 
F(k;t) 



dS(k) S(k) 



dk 

5 dS(k) 
k dk 



D a k 2 S(0)nc(k) 
6 



k 
2 

d 2 k F(k;t) 



dS(k) 
dk 



-i 2 



S(k) S(k) 2 dk 

2 fdS(k)\ 2 d 2 S{k)~ 
5(fc) \ dk J + dk 2 

4D k 2D k 2 dS(k) 



D a k 



S(k) S(k) 2 dk 



S(k) 
d k F(k;t) 



F(k;t) 



Dok 2 
' S(k) 



d 2 F{k;t) 



~d 2 S{q) 




5(0) ~ 




dq 2 


q=0 


k 2 


I 



dt'M ilr (k;t-t')d t 'F(k;t') 



5(0) 



dt'M m (k-t-t')d t id k F(k\t') 
dt'M^{k-t-t')dt>d 2 k F{k-t'), 



and 



M*(k;t) = nD 



dki 



[^(k^k-ki)] 5 



(2tt) 3 1 
Fdk-kil;^^!^), 



and 

M*(fc;t) = 



(22) 



(23) 



j^ } v k (k 1 ,p)F(p;t) X {1) (k 1 ;t) 



(2tt) 3 

^(ki.p) _ 2(k-ki) 

k k 2 k x 



«k(ki,p) 



k • ki dc(ki) 
k dki 



uDq f dki 



«k(ki,p)F(p;t)x (0) (*i;<) 



2 J (2tt) 3 
c(fci) _ ^k(ki,p) k ■ ki d 2 c(fci) 
jfe "~ k 2 2k dk\ 

2(k • ki) 2 - 2fc 2 k • ki - k 2 k\ dc(fei) 
fc 3 fci dfci 



In Eq. (231 p = k - k x and p = |k - kj|. The initial 



condition is given by 

x (2) (fc;0) = ^ 



d 2 5(q) 



dq 2 



9=0 



(24) 



1 5(fc)5(0) dS(k) 1 d 2 5(fc) 
+ 3^fc ^ + 6 5m ° ) ^fc 2 - 

To solve Eq. [20] we also need the equation of motion 



(17 1 for d k F(k;t) and the equation of motion for the 
second partial derivative of the intermediate scattering 
function with respect to the wave- vector, d 2 F(k\ t). The 
latter equation can be obtained from the mode-coupling 



+ / dt' M' m '(k;t - t')d r d 2 k F{k-t') 



dt'Ml 2 (k;t-t')d t id k F(k; t') 



+ / dt'M« 2 (k;t-t')d t ,F(k;t') 
Jo 

= 0, 
where 

M* 2 (k;t) = 

9 n f dkl n, J < u P^Pl!^)! 
2nDo J (2^ Wk(kl ' P) \ C{P) + l^p-^p-) 



(25) 
(26) 



xF(k i; t)F(p; t) 
dki 



-nDo 



( 3 [ Wk (k!, p)] 2 ^9 p F(p; * - t')F(kv, t) 



and 

rfe2 



(fc;i) = 



(27) 



nD 
+nD 



dki 
f^rj 3 
dki 



F(k i; t)F(p;t) 



c(p) 



[k ■ p] 2 dc(p) 
fc 2 p dp 



-i 2 



(2tt): 



-u k (k 1)P )F(fci;t)F(p;t) 



3k ■ p dc(p) _ [k ■ p] 3 rfc(p) + [k ■ p] 3 d 2 c{p) 
kp dp k 3 p 3 dp k 3 p 2 dp 2 



-2nD c 



dki 



«k(ki,p) 



c(p) 



[k-p] 2 dc(p) 



xF(kf,t)^-d p F(p; t) 



nDn 



dki 



2 7 (2tt) 3 
x9 pJ F(p; t)F(h;t) 

nD f dki 

2 7 (27T) 3 



Mki, P )] 2 

w k (ki,p)] 2 



/c 2 p dp 
1 [k-p] 2 " 



p k 2 p 3 



k P 

kp 



F(fci;t)a 2 F(p;i) 



In Eqs. ([26J27) p = k - k x and p = |k - ki| 



IV. ISOTROPIC APPROXIMATION 

Along with the expansion of the full equation of mo- 
tion Q , we also examined an expansion of an isotropic 
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approximation to Eq. Q. The isotropic approximation 
has the advantage of being slightly easier computation- 
ally, and it allows for calculation of the susceptibility at 
any q. 

The isotropic approximation assumes that the suscep- 
tibility Xq(k; t) is independent of the angle between q and 
k, Xq(k; t) ~ x!j SO (^; £)• To derive an equation of motion 
for Xg°(^j t) one could start by substituting the isotropic 
approximation into the full equation of motion and then 
average the resulting equation over the angle between q 
and k. We propose a slight modification of this proce- 
dure that results in an equation that is somewhat easier 
comput ationally : 



dt* i q so (k ] t) + 



S(k) 



xL so (M) 



+ / dt'M irr (M-0^'Xa so (M') 



+ / dt'M™ {k;t-t')dt>F(k;q;t') 



= nD k 2 S(0)c(k)F(k;t) 



-5(0) / dt'M' m '(k;t~t')d t 'F(k;t'). 
'o 



(28) 



where 

M™{k;t) 



nD, 



dkx 



X^fciJtjFflk-kxl;*) 



(27r) 3/V<J 

xi> k (ki,k - ki)« k (ki,k - ki; q) 



and 



F(k;q;t) 



dq 

4tt 



F(|k + q|;t), 



«k(ki,k- ki;g) 



dq fcu k +q(ki + q, k - ki) 
4tt |k + q| 



(29) 



(30) 



(31) 



Finally, the initial condition for X« so (k;i) is 



^(k;t = 0) = S(k)S(q) 



dq 

47T 



5(|k+q|). (32) 



Note that in Eq. (28 1 we took the source term in the 



q — > limit. The q dependence of the source term has 
very little effect on the size of the correlation length (see 
discussion in Sec. VI I. Taking the source term in the 
q — > limit makes the numerical calculation of Xg°(ki t) 
somewhat easier. 

To get the characteristic length we expanded Xg S °(^i t) 
in powers of q, 



xf(M) - x (0) (M) + 



dq 2 



J q=0 



x (0) (M) + <z 2 x iso(2) (M) 



(33) 



The zeroth order coefficient, X^°\k',t), is the same 
as the one obtained from the expansion of the com- 
plete equation of motion. The equation of motion for 
X lso( - 2 )(/c; t) can be readily obtained from Eq. (28 1: 



d t x iso{2) (k;t) + ^x iso{2 Hk;t) 



+ / dt'M 1TI (k;t-t')d t/ x lso{ ) (k;t) 

JO 

+ / dt 1 'M£°(k;t - t')d t >F(k;t') 



dt'M$(k;t-t')d t >d%.F(k;t') 
di?M$(k;t-t')d t ,d k F(k;t') 



6 Jo 



3& Jo 



dt' 'Mf°{k;t - t')dt'F(k-t') 



where 
M™°{k;t) 



(34) 



(35) 



M x (fc;t) is defined in Eq. (13 1, and 



MT(M) = 

^ J ^M^,p)F( P -,t)x^(ki;t) 

" c(fci) _ T; k (ki,p) k ■ ki d 2 c(ki) 

k k 2 2k dkj 

2(k • ki) 2 - 2fc 2 k ■ ki - k 2 k\ dc{k t ) 



(36) 



In Eq. (361 p = k-ki and p = |k — k x | 



The initial condition to Eq. ( 34 1 is given by 
S 2 (k) 



X iso{2) (k;0) = 



d 2 S(q) 



dq 2 



q=0 



(37) 



S(0)S(k) dS(k) S(0)S(k) d 2 S(k) 
3k dk~ + 6 dk 2 ' 



V. NUMERICAL EVALUATION OF x ( '°(M) 

We numerically calculated the k and t dependence 
of X^°Hk;t), X^{k;t) and X^ 2 \k;t) using a previously 
developed algorithm that was designed to solve mode- 
coupling like equations |15[ 116} I17j . The only input in 
this calculation is the static structure factor S(k), which 
we calculated for the hard sphere interaction potential 
using the Percus-Yevick approximation. We report our 
results in terms of the relative distance from the ergodic- 
ity breaking transition predicted by mode-coupling the- 
ory, e = (<f> c — 4>)/4> c . Here <fi IS the volume fraction, 
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e dependence, i.e. tmix ~ e 2A6 . Shown in the Fig.^ are 
the ratios imax/T Q , and it can be seen that these ratios 
are constant close to the mode-coupling transition. Thus, 
in the e — » limit we see that imax = 1-4t q , imax = 2.8t q , 
and tmL = 1.4r a . Notice that the peak of X^C^max!*) 
and x^ 2 H^i 

Finally, in Fig. [3] we compare the values of 



( (n) (fcmax;i) att 



V^rnax, ^max 
3/2 



max j 1 

max, £) occur at the same time, 
in Fig. [3] we compare the 
We find that | X (0) ( 

grow as e -1 whereas |x (2) (^max; fmaxj 
As we discuss in the next section, this 



(n) 

max- 



(0) 



and 

(2) 



FIG. 1: Time and wave- vector dependence of q = value 
of three-point susceptibility, Xq=o(k;t) = \ (k;t) for the 
reduced distance from the mode-coupling transition, e = 
(4> c — 4>)/4>c = 0.05 (upper panel) and e = 10~ 4 (lower panel). 
Contours correspond to x (k;t) = 4 n where n is an integer, 
starting from n = —6. The arrow marks the position of the 
first peak of the static structure factor. 



<f> = nira 3 /6, where a is the hard sphere diameter, and <p c 
is the volume fraction at the mode-coupling transition. 
We used 300 equally spaced wave-vectors with spacing 
S = 0.2, between k — 0.1 and fc max = 59.9, and this 
discretization resulted in a mode-coupling transition at 
C = 0.515866763. 

Shown in Fig. [T| are contour plots of x^ Q \k;t) as a 
function of wave-vector k and time t for e = 0.05 and 
e = 10~ 4 . The former value of e is the smallest relative 
distance from an avoided mode coupling transition in the 
Kob-Andersen binary mixture at which mode-coupling 
theory agrees with computer simulations |17j . As we 
mentioned earlier, x^°\k;t) is proportional to the three- 
point susceptibility Xn(k;t) calculated in Ref.[18 j , which 
is a mode-coupling approximation for the density deriva- 
tive of the intermediate scattering function. Thus, all 
results derived in Ref.[T5] for for Xn(k)t) also apply to 



X 

Ix 

grows as e 

disparate behavior of X (fenaxi t) is important for the 
existence of a diverging characteristic length. 

We should note at this point that the e dependence 
of x (n) (k;t) can be deduced from scaling predictions de- 
scribed in Rcf. |13j . and the numerical results presented 
here fully agree with the these predictions. 



VI. DIVERGING CHARACTERISTIC LENGTH 

To obtain a growing characteristic length scale as 
the mode-coupling transition is approached, we need 
\x^ n \k; t)\ for some n > to grow faster than |x'°H^ *)l 
at a fixed time t. Then a diverging length can be calcu- 
lated as |x (n) lfc;*)/x (0) (fc;*)l (1/n) - 

From Fig. K3nt is clear that the linear term, X *)> 
does not result in a growing length scale. On the other 
hand, the absolute value of the isotropic second order 
term, \x^ (k;t)\, grows faster than x (k;t) and thus 
we can define a diverging characteristic length £(fc; t), 



£(M) 



x (2) (fc;*) 
x (0) (M)' 



(38) 



(0) 



where the negative sign comes from the observation that 
X^ 2 \k;t) is of opposite sign of x^°\k;t) around r a and 
close to the transition. Note that for large times t, Eq.[38| 
involves a division of a small number by another small 
number. Because of numerical issues present in the al- 
gorithm to calculate x^ n \k\t), we only show results if 
7k; t).~In particularTtnere is a well denned maximum * (n) ^ l ) ^ 10 ~ 3 ' and therefore we, unfortunately, can 



X 

in x (ty t) at a well denned wave- vector and at a char- 
acteristic time. Also, all the scaling laws observed for 
Xn(k;t) apply to x (k',t) (we show some of these scal- 
ing laws below). The characteristic wave- vector, fc maX) 
is nearly constant as the mode coupling transition is ap- 
proached and fcmax ~ 7.1 close to the transition. 

While the characteristic wave- vector is nearly constant 
close to the transition, the characteristic time grows 
rapidly as the mode-coupling transition is approached 
and diverges at the transition. In Fig. [2] we examine the 
time at which x^ (tai! t) is a maximum, tmax, as a func- 
tion of e for the characteristic wave- vector fc max = 7.1. 
We compare t 



not comment at the asymptotic t — > oo limit of the char- 
acteristic length. 

In Fig. [4] we examine ^(fc max ; t q ), i.e. the characteristic 
length at k = fc max and at the a relaxation time. The 
length £,(k max ;T a ) grows as e -1 / 4 and it reaches only 15 
particle diameters at e = 10~ 6 . Thus the characteristic 
length is not very large even very close to the transition. 
For e = 0.05, we find that £(fc max ;T Q ) is only about one 
particle diameter. Note that Eq. (38 1 defines a length 



(n) 



max with the a relaxation time r a , for which 



scale for every wave- vector k and at all times t, and we 
examine the time and wave- vector dependence of £(k;t) 
below. 

We determined that setting the initial condition for 



we use the standard definition F (fc ma x; T a ) 



We 



X 



(2) 



(n) 

find that i ma x is slightly larger than r Q , but has the same 



(k;t = 0) to zero and/or taking S^(k;t) = had 
very little effect on the size of the correlation length close 
to the mode-coupling transition. While including these 
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FIG. 2: Upper Panel: the a relaxation time, r a (filled circles) 
and the peak positions of X 



((f> — (f>c)/4>c- Tmij 



T-max, as a function of 

-triangles; Tmax-diamonds; r4 2 2 x -open 

(n) 



squares. Lower Panel: the ratio Tm&x/Ta as a function of e 




FIG. 3: The peak height of x^ n \k;t) as a function of the dis- 
tance from the mode coupling transition e. X (^max! tmax) ~ 
triangles; x (1) ( fc ma X ; imax) - circles; x (2) (fcma X ; tmax) - squares 



terms is in principle straightforward, dropping them sig- 
nificantly simplifies the numerical calculation. 

Fourier transforms of four-point correlation functions, 
i.e. four-point dynamic structure factors, are often mon- 
itored in simulations and used to investigate properties 
of dynamic heterogeneities. Since the q = value of a 
four-point structure factor should be proportional to the 
characteristic volume in which correlated motion takes 
place, an increase of the q = value (i.e. of the height 
of four-point structure factor) is often given as evidence 
of an increase in a dynamic correlation length. 

Similarly, for the problem considered here, the value of 
x(°)(fc; t) could used as an indicator of the size of a char- 
acteristic dynamic range of the response. However, the 
spatial extent of dynamic response is best measured by 
examining the long-range spatial decay of a direct space 



FIG. 4: The characteristic dynamic length £(fc ma x) T a ) as a 
function of the distance from the transition e. 



susceptibility or, alternatively, by examining the small-g 
behavior of the susceptibility Xq(k;£). This distinction 
is significant in view of the very strong wave-vector and 
time dependence of x (k;t). In particular, if the char- 
acteristic length were a monotonic function of x^ -* (fc; t), 
then Fig. [I] would be leading to the unfortunate conclu- 
sion that £(fc; i) is a very strong function of k. The length 
would then have a rather limited appeal. In the following 
paragraph we show that this is not the case. 

In Fig.[5]we compare the k dependence of £(£;; t) (right 
figure) and x °\k;t) (left figure) for three characteris- 
tic times: (1) early (3 (dotted line), late-/3 (dashed line), 
and at the a relaxation time (solid line). For reference, 
F(k max ; t) is shown in the insert to Fig. [5] with the three 
characteristic times shown as vertical lines in the figure. 
There is a very strong dependence of x {k; t) on k, but 
£(/c; r a ) is nearly constant at each time. Therefore, even 
though there is a strong k dependence of the three-point 
susceptibility, there is a well defined characteristic dy- 
namic length £,(k;t) that is independent of k and only 
depends on the time t. This suggests that we could drop 
the explicit k dependence of £(fc; t) and introduce a sim- 
plified notation £(t). 

Next, we investigate the time dependence of the char- 
acteristic length. Shown in Fig. [H] is x^H^max', t) (lower 
curve-right axis), Ix^f^max; £)| (middle curve-right axis) , 
and £(£) (upper curve-left axis) as a function of time for 
e = 10~ 6 . The correlation length £(t) is close to one for 
t = 0, begins to grow during /3 relaxation and reaches a 
plateau at a time corresponding to the late /?-early a re- 
laxation. During the a relaxation, £(i) is approximately 
constant. Note that £(t) has a very different time de- 
pendence than x^H^max; t). Therefore, the length scale 
associated with dynamic heterogeneities are not a max- 
imum when x (0) (k; t) is a maximum, but rather reaches 
a constant value for times less than this characteristic 
time. 

Scaling relations for different time regimes can be de- 
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FIG. 5: Left panel: the wave-vector dependence of the q = 
value of the three-point susceptibility, Xq=o(k;t) = x [k;t), 
at a time corresponding the the early relaxation regime 
(dotted line), the late j3 regime (dashed line), and the a relax- 
ation time (solid line). Right panel: the wave- vector depen- 
dence of the characteristic dynamic length £(fc; t) for the same 
times as in the left panel. The inset is the self-intermediate 
scattering function F(k; t) and the three vertical lines corre- 
spond to the three times in left and right panels. 





Mill | Hill | Ilium 


i i i 




| lull | limit 


10 1 


: e=l(T 6 


a/2 '/ 














/y' b \ 








2a 












M 








10° 


■ ffsr 


i i i 







r 



10 

10 6 1 
10 5 J 

io 4 C 
10 3 ?5 

110% 
110° 



10 



10" 



FIG. 7: The characteristic dynamic length ^(fc max ;i) as a 
function of t/r a for e = 0.05, 10~ 4 and 10~ 6 . The dashed 
lines is the scaling law £(t) ~ t a ^ 2 valid in the j3 relaxation 
regime. 
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FIG. 8: The characteristic dynamic length £ (r a ) calculated at 
the a relaxation time as a function of the a relaxation time. 



FIG. 6: The time dependence of the characteristic dynamic 
length £(t) (left solid line and left axis), the susceptibility 
X^°' (^max! t) (right solid line and right axis), and the second 
order coefficient X (^maxji) (middle, heavy solid line and 
right axis), showed on a log-log scale. The dashed lines show 
the scaling laws in the f3 relaxation regime. The vertical lines 
crossing £(t) correspond to the three times shown in the inset 
to Fig. [5] 



rived from the predictions of the mode coupling theory 
[5JH3]. Specifically, in the early (3 regime X (^i*) ~ 
and in the late /3 regime x^H^i ^) ~ t b where a = 0.312 
and b — 0.583 for our system. The power law growth 

of X^(*) an d X^ft) are a -^ so snown in Fig- ® Dur- 
ing the early (3 relaxation regime, X (t) ~ t a while 
X^Ht) ~ t 2a , which gives rise to the W 2 growth of the 



correlation length in the early (3 relaxation regime. How- 
ever, during late (3 relaxation, x' 2 '(*) an d X^°\t) both 
grow as t , thus there is no growing length scale. The 
vertical lines in the figure denote the same times as the 
vertical lines in the inset to Fig. [5] 

In Fig. [7] we show £(fc max ;£) as a function of t/r a for 
e = 0.05, 10" 4 , and 10~ 6 . For e = 10~ 4 we observe the 
t a / 2 scaling for only a very narrow range of time, and 
we do not observe the W 2 scaling for any time range at 
e = 0.05, which suggests that it might be very difficult 
to see this scaling in simulations. 

Finally, we note that since r a ~ e - 2 - 46 an d £ ~ e -0 - 25 , 
then £ ~ r°' 102 , Fig. [8j As a result, a modest increase 
in the correlation length is accompanied by a very large 
increase of the relaxation time. 
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VII. NUMERICAL EVALUATION OF x£° (M) 
AND ASSOCIATED CHARACTERISTIC 
LENGTH 

The isotropic approximation neglects the dependence 
of the three-point susceptibility on the angle between k 
and q. Thus, in the resulting equation of motion for 
X 1 9 so (fc; t) q is just a parameter, and the equation of mo- 
tion can be solved separately for any value of q. As a re- 
sult, the full q dependence of Xg SO (fcj t) can be calculated. 
On the other hand, the isotropic approximation preserves 
the essential terms in the equation of motion which lead 
to the divergence of the q — > limit of Xq°(k', t) and of 
the characteristic length. In this section we examine the 
isotropic approximation and compare this approximation 
to the expansion terms given above. Since the equations 
of motion are similar and the terms that cause the diver- 
gence are identical, many of the results of Sec. |Vj carry 
over to the isotropic appr oximation . Notably, as we al- 

X^(k;t) is identical in both 



IV 



ready noted in Sec. 
cases. 

Since we can calculate the whole q dependence in the 
isotropic approximation, we can determine the character- 
istic length using two different methods. We can ei- 
ther evaluate Xq (k; t) and then fit Xq°(k; t)/x {0) {k\ t) to 
1 — (£ lso (&:; t)q) 2 for small q or we can determine £ lso (fc; t) 
from ^-x lso(2) (k; t) /x (0) (k; t) . Both methods result in 
the same length. 

It can be showed that within the isotropic approxi- 
mation the characteristic length is almost /c-independent 
(and thus we will denote it by £ lso (t))- In addition, the 
time dependence of the length is very similar to what was 
obtained from the full equations of motion in Sec. [V] 

In Fig. |] we compare the magnitude of the charac- 
teristic length obtained from the isotropic approxima- 
tion, £ 1so (tq,), with that following from the full equations 
of motion, £(t q ). As we anticipated in the first para- 
graph of this section, the isotropic approximation gives 
a length which diverges as e -1 / 4 . However, the isotropic 
approximation underestimates the characteristic length; 
for small e the length resulting from the isotropic approx- 
imation is approximately 36% smaller than the length re- 
sulting from the expansion of the complete equation Q . 

There has been some discussion in the literature as to 
what scaling function should be used to determine £(t). 
According to the scaling relation presented in Ref. |13j, 
in the /3 and a regimes the divergent part of Xq(k; t) is a 
function of a scaling variable q£(t) only for small q close 
to the transition, Xq(k;i) = <^8,a(<z£(*|9,a)» k). We use 
the isotropic approximation to examine some properties 
of scaling function Xp^ a close to the mode coupling tran- 
sition, in the f3 and a regimes. 

For times t in the vicinity of the (3 relaxation time Tp , 
the scaling function Xp(q£(tp), k) is predicted to have the 
Ornstcin-Zcrnicke behavior, namely Xp{Q£{tp), k) should 
scale as q~ 2 for large q [T3]. To check this prediction we 
first need to define the (3 relaxation time. We define 773 as 




FIG. 9: The characteristic dynamic length £ lso (fc; t) calculated 
using the isotropic approximation (squares) and without the 
isotropic approximation (circles). 




FIG. 10: The isotropic approximation for the dynamic sus- 
ceptibility Xg°(fcmax; 7/3 ) as a function of q£(Tp) for e < 10 -3 . 
Only data for q in the scaling regime are included. The solid 
line is the Ornstein-Zernicke function 1/[1 + (£q) 2 ]. 



the inflection point of F(t) versus ln(i). We verified that 
this definition agrees with the MCT scaling Tp ~ e~ 1//2a . 
This time Tp is only well defined for e < 10~ 3 . Shown 
in Fig. 10 is x™(^max; Tp) /x (0) (W; Tp) as a function of 
q£,(T a ) and the Ornstein-Zernicke function 1/[1 + (£<?) 2 ], 
which provides a good fit for small q during the [3 relax- 
ation time and demonstrates the q~ 2 scaling for large q. 



For times t comparable to the a relaxation time r Q , 
the inhomogeneous mode-coupling theory |13j predicts 
a q~ 4 behavior of the scaling function X a (q£(t a ) , k)&t 



large q. 

max * 



We test this prediction in Fig. [TT] we show 
q vmax,T Q )/x (0) (fcmax;T o ) as a function of q£ (r Q ) along 
with two functions commonly used to find £(t) in simu- 
lations, and a function suggested by the inhomogeneous 
mode-coupling theory. The functions 1 — (£q) 2 (dotted 
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FIG. 11: The isotropic approximation for the dynamic sus- 
ceptibility X9°(fcmax; T a ) as a function of q£,(r a ) for e < 10~ 3 . 
Only data for q in the scaling regime are included. The dotted 
line is 1 — (q$,) 2 , the dashed line is l/[l+(g£) 2 ]. The solid line is 
a fit to the data to a function of the form 1/[1 + {qO 2 + a (?£) 4 ] 
where a = 0.45. The q~ 4 scaling for large q is evident. Inset: 
the dynamic susceptibility Xq °(&max; r a ) for e = 0.05 showing 
all the data including q values beyond the scaling regime. The 
inset shows that the q~ 4 scaling (solid line) is not apparent 
for this e. 

line) and the Ornstein-Zernicke function, 1/(1 + [£(/] 2 ), 
(dashed line) are good fits only to a very narrow q range, 
with the Ornstein-Zernicke function being a better fit for 
a larger range of q values. On the other hand, the func- 
tion l/[l + (^q) 2 +a(^q) 4 ] where a — 0.45 (solid line), pro- 
vides a good fit over a large q range and thus it confirms 
the q~ A scaling predicted by the inhomogeneous mode- 
coupling theory for the a relaxation time scale. Note that 
the q~ 4 scaling is not evident for e — 0.05 (inset), which 
suggests that this scaling might be difficult to observe in 
simulations. 



VIII. CONCLUSIONS 

We used inhomogeneous mode-coupling theory to nu- 
merically investigate the dynamic susceptibility Xq(k; t) 
at small q and determined time, fc, and distance from 
the transition dependence of the diverging characteristic 
length scale. We confirmed scaling predictions presented 
in Ref. [T3] and obtained a couple of new interesting 
results. Because of numerical issues, we were not able 
to calculate the asymptotic long time behavior of the 
diverging characteristic length scale. This would be an 
interesting topic that we leave for later analysis. It most 
likely requires an analytical argument that goes beyond 
the scaling analysis presented in Ref. [13] . 

The most important result of our numerical investi- 
gation is that the diverging characteristic length is very 
weakly k dependent. This makes it a well defined quan- 
tity. We speculate that the k independence of the char- 



acteristic length should carry over to the dynamic cor- 
relation length defined in terms of a four-point struc- 
ture factor. Moreover, it should explain why a variety 
of slightly different four point functions (e.g. defined 
in terms of overlap functions [HE [2U] or in terms of 
self-intermediate scattering functions [6j|7j|2T]) result in 
comparable dynamic correlation lengths. 

The second important result, which cannot be obtained 
from scaling considerations alone, is the magnitude of the 
characteristic length. On general grounds we expect this 
length to be comparable to dynamic correlation lengths 
defined through four-point structure factors. Thus, it 
is satisfying that the magnitude of the length is indeed 
comparable (albeit somewhat smaller) to what's found in 
simulations. 

We would like to point out that, although various sim- 
ulations found comparable values of the dynamic cor- 
relation length, there are a few important unresolved 
differences between results obtained by different groups 
that preclude declaring that the characteristic length dis- 
cussed in this work is essentially the same as the dynamic 
correlation length measured in simulations. 

First, while the characteristic length defined through 
the three-point susceptibility is a monotonic function of 
time (at least as long as our numerical routines are re- 
liable), the simulational results very. Lacevic et al. [5] 
found that the dynamic correlation length roughly fol- 
lowed the overall magnitude of the four-point correlation 
function and decayed to zero at later times. In contrast, 
Toninelli et al. [7] found that the dynamic correlation 
length continued to grow at later times. While slightly 
different fitting procedures were used in these two works, 
it is difficult to pinpoint the exact source of two strikingly 
different results. 

Second, within the mode-coupling approximation, the 
characteristic length defined through the three-point sus- 
ceptibility diverges as e -1 / 4 upon approaching the ergod- 
icity breaking transition predicted by the mode-coupling 
theory. We feel that the relevance of this result to sim- 
ulations (and experiments) in which the mode-coupling 
transtion is avoided still needs to be fully established. We 
speculate that it is possible that in computer simulations 
a vestige of a power law divergence of the dynamic cor- 
relation length could be seen just as one can observe in 
simulations power law dependencies of various transport 
coefficients upon approaching a mode-coupling crossover 
[T7] . Indeed, various groups have already claimed power 
law dependencies of their dynamic correlation lengths 
upon approaching the mode-coupling crossover (see, e.g. 
[SI OH HO HI])- However, there seems to be some 
disagreement regarding the value of the scaling exponent 
and only one work, |19j . results in a value agreeing with 
the prediction of the inhomogeneous mode-coupling the- 
ory. Upon closer examination of the fitting procedure de- 
scribed in Ref. [23 and re-examining our own simulation 
data we concluded that virtually all systems studied in 
simulations were not large enough to obtain the dynamic 
correlation length in a range allowing for an unambiguous 
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